#include <opencv2/gpu/gpu.hpp>
#include <highgui.h>
#include <iostream>
#include <cmath>
#include <cv.h>
#include "mi.hpp"

# define DIM 256

using namespace cv;
using namespace std;

double calculate_Mutual_Information(Mat & A, Mat & B)
{

   double entropyA,
		  entropyB,
		  jointEntropy;
   Mat	  jointHistogram,
		  jointProbabilityDistribution;


	entropyA      =  calculate_Entropy(A);
	entropyB      =  calculate_Entropy(B);
	jointEntropy  =  calculate_Joint_Entropy(A,B);

	//cout<<"Joint Entropy:  "<<jointEntropy<<endl;
	//cout<<"A's Entropy:    "<<entropyA<<endl;
	//cout<<"B's Entropy:    "<<entropyB<<endl;

	return entropyA + entropyB - jointEntropy ;

}

Mat calculate_Probability_Mass_Function(Mat & image)
{

    float range[] = { 0, 256 } ;
	const float* histRange = { range };
    const int histSize = 256;
    bool uniform = true; bool accumulate = false;

	Mat  normHist,
	     hist;

	normHist.create(hist.rows, hist.cols,CV_64F);
	calcHist( &image, 1, 0, Mat(), hist, 1, &histSize, &histRange, uniform, accumulate );
	normHist.create(hist.rows, hist.cols,CV_64F);
	normHist =  hist / (double)(image.cols * image.rows);

   return normHist;
}

double calculate_Entropy(Mat & image)
{
	Mat     normHist;
	double  probability = 0,
		    entropy     = 0;

	normHist = calculate_Probability_Mass_Function(image);

	for(int i = 0;  i < normHist.rows; ++i){

		for(int j = 0; j < normHist.cols; ++j){

			probability = normHist.at<float>(i,j);

			if(probability > 0.0)
			{
				entropy += probability * log2(1/probability);
			}
		}
	}
	return entropy;
}

Mat create_Joint_Histogram(Mat & A, Mat & B)
{

	Mat jointHistogram;

	jointHistogram = Mat::zeros(DIM, DIM, CV_32FC1);

	for(int i = 0;  i < A.rows; ++i){

		for(int j = 0; j < A.cols; ++j){

			jointHistogram.at<float>(A.at<uchar>(i,j),B.at<uchar>(i,j))++;
		}
	}
	imshow( "jointHistogram", jointHistogram);

	return jointHistogram;

}

Mat calculate_Joint_Probability_Distribution( Mat & A,  Mat & B)
{
    Mat	jointHistogram = create_Joint_Histogram(A,B);
    Mat jointProbabilityDistribution =

    		 jointHistogram / (A.cols * A.rows);

 return jointProbabilityDistribution;
}

double calculate_Joint_Entropy( Mat & A,  Mat & B)
{
	double  jointEntropy = 0,
			probability  = 0;
    Mat     jointProbabilityDistribution;

	jointProbabilityDistribution =
		calculate_Joint_Probability_Distribution(A,B);

	for(int i = 0;  i < DIM; ++i){

		for(int j = 0; j < DIM; ++j){

			probability = jointProbabilityDistribution.at<float>(i,j);

			if(probability != 0){
				jointEntropy += probability * log2(1/probability);
			}
		}
	}
	return jointEntropy;

}









